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ANALYSIS OF A FORTRAN PROGRAM 
FOR COMPUTING ELECTRIC FIELD DISTRIBUTIONS IN HETEROGENEOUS 
PENETRABLE NONMAGNETIC BODIES OF ARBITRARY SHAPE 
THROUGH APPLICATION OF TENSOR GREEN'S FUNCTIONS 


PURPOSE 


TGR is a Fortran IV program written for the IBM 360/65 System to 
compute the electric field distribution, power density distribution, 
and total absorbed power in a bounded, penetrable, heterogeneous, non- 
Magnetic body subject to an incident electric field. Results are 
obtained through the discretization and simplification of the integral 
equation of electromagnetic scattering derived from Maxwell's equations 
by applying simplifying assumptions appropriate for the method of 
retarded potentials. 


In the body of the text we give a terse description of the 
integral equation that we use. Appendix A contains the mathematical 
preliminaries needed to understand the complete derivation, which is 
given in Appendix B. The discretization of the integra! equation 
and its representation as a matrix equation is explained in Appendix C. 
The reduction by the use of symmetry groups of the size of the matrix 
used to define the approximation of the integral equation is explained 
in Appendix D. Without using the latter technique the cost of running 
the program would be greater by from one-hundred to seven-hundred 
percent than it is at the present time when the scattering obstacle 
possesses a lesser or a greater amount of symmetry. 


With this program the authors were able to observe resonance 
effects and predict hot spots in a biological body subjected to 
microwave radiation. In this preliminary report a comparison is 
made between results produced by the method of this report and the 
classical Mie solution results in the case where the impinging field 
is a plane wave and the scattering obstacle is a homogeneous, non- 
biological ball with a relative dielectric constant close to unity 
and a conductivity that is small. Good agreement was found between the 
computations produced by the two methods; the results of two numerical 
experiments are reported in Appendix E. The error analysis for biological 
bodies is being carried out, and the results will appear in a future 
report. The listing of the computer program along with the computer 
output for one of the sample problems is given in Appendix F. 


MATHEMATICAL DESCRIPTION 


A biological body is subjected to a time-harmonic electromagnetic 
wave traveling in the direction of the positive z-axis. This incident 
electric vector is normalized so that its amplitude is 1 volt per meter. 


The incident electric vector may have the form 


t x4 cos(kpz)s (1) 

EY = ¥ sin(kyz), (2) | 
or 

EY = ¥ (cos(kgz) - isin(kyz)), (3) 


where { is a unit vector in the direction of the positive x-axis ’ 
and kg = wlgeg is the free space wave number. 


We use the MKS (meter-kilogram-second) system of units. The magnetic 
field intensity H which is measured in ampere-turns per meter is given 


by 

A= (1/iwny)(2/azEi , (4) 
or | 

Ho = - vVleq/ug)€,exp(-ikgz). (5) 
The Poynting vector (Stratton [6 , p. 132]) is given by 

S$ = ExA, (6) 
and has the units of watts per square meter. In the MKS units 

f° 8.854x10"!* farads per meter, - (7) 

ug = 1-257x107° henrys per meter, (8) 
and 

VeqTig = 2-654x107> mhos - (9) 
Thus, the power density for a 1 volt per meter electric field is 

(sj = 2.654x10°> watts per square meter . (10) , 


According to Johnson et al. [1, p. 14], power density is normally 

measured in milliwatts per square centimeter. If the amplitude of . 
the incident electric field were multiplied by 61.38, then the 

power density would be 1 milliwatt per square centimeter. 


a 


To get the power that would be absorbed from a 1 milliwatt-per-square- 
centimeter beam, we simply multiply the normal output of this program 
by 3767., remembering that the wave impedance of free space is 376.7 
ohms according to Stratton (6, p. 601] . 


The authors define the tensor Green's function 


(-ik =F) 
G(r.r') = (17 (4n)) {1 + (ME rarg grag) | Pol (11) 


[r= | 


and show the connection between the Maxwell equations, the calculations 
of Van Bladel {7, p. 563] and Livesay and Chen [4, p. 1273], and obtain by 
neglecting surface electric current and surface magnetic charge the inte- 
gral equation 


[ + Jer) : Pv [c(t )E(r* G(r ar )av(r * € tr). (12) 


where the PV indicates that the principal value of the integral is taken, and 
t(r) = o(r) + iw(e(r) - Eq): (13) 


Equation (12) is solved by subdividing the body into subregions in which the 
electric vector is assumed to be constant as a function of space coordinates. Thus, we 
assume, at least, that the diameter of a subregion does not exceed one-fourth of a 
wavelength. Since the velocity of light is about 30 billion centimeters 

per second, this means that at 1 megahertz the diameter may be only 7.5 
centimeters and at 1 gigahertz the diameter may be only .0075 centimeters. 
Thus, without modification, it would be impossible on the present computer 

system to determine accurately the effect of 1 gigahertz radiation on a 
180-centimeter man, because of the memory that would be needed to store 
bai about all of the cells that would be involved in the approxima- 

tion. 


PROGRAM DESCRIPTION 


The program is adapted from a program written by Livesay and Chen [4]. 
We have adopted the same input and output format adding only the capability 
of printing the total absorbed power. This program is faster than Chen's 
program and requires less storage in view of the fact that the matrix 
inversion routine stores and operates only on about half of the symmetric 
matrix defined by the discretization of (12). 


Following Chen we divide space into octants. The first octant is that 
region of x-y-z space where x,y and z are all positive. Part or all 
of the scattering obstacle is assumed to occupy octant one. If it occupies 

any other octant it is symmetric with respect to that part that is in octant 
one. The use of symmetry is explained in Appendix D. 


We now explain the input data. Table 1 lists the files for this input 
data. Table 2 is a listing of the formats for these files. 


TABLE 1. INPUT DATA DESCRIPTION 


Data File Contents 


N DIV 


al COMP: Q(J), J=1,...,8; FMEG, FIELD 
NX, NY, NZ 
N 


ao fF WN 


{AMX, AMY, AMZ, RLEPI, SIGI, DXCM}, 


j 
N cards 


{AMX, AMY, AMZ, RLEPI, SIGI, DXCM} 


qf Here, NDIV contains the number of subdivisions of each side of a subvolume. 

Next, COMP is a code name for the components of the induced electric field; 

1 it contains one of the members of the set of alpha strings denoted by 

: {X-ONLY, XANDY., X, Y, Z.}. If an octant J is used, then the number of that 

| octant is placed in the field for Q(J) and otherwise a blank is placed in the 
field; octant 1 is always used. The variable FMEG denotes the frequency in 
megahertz. The variable FIELD specifies the character of the input field and 

is an a-string belonging to the set {EXPKZ, COSKZ, SINKZ, TOTAL, EINCZ} for 

i | respectively fields of the form 


EXPKZ = 7 


— 


exp(-ik.z), (if we use only octants 1 - 4) 


i 
i COSKZ 


=_ cos(koz), 
SINKZ = 7 sin(k.z), 
: . 
i TOTAL = 7 exp(-ik9z). 


EINCZ 


The variables NX, NY, and NZ approximate the number of cells in the X, Y, 
and Z directions; this information is not used. The variable N denotes the 
number of cells. The variables AMX, AMY, and AMZ correspond to the maximum 
values of X, Y, and Z coordinates of a cell in the first octant. The 
variables RLEPI and SIGI denote the relative dielectric constant and 
conductivity in mhos per meter of the cell in question. The variable DXCM 
is the length of the edge of the cubical cell in centimeters. The values of 
RLEPI and SIGI may be determined from [1] and [2]. 


I. The output is clearly labeled by the program and is self-explanatory. 
The reader should be able to see how the program is used by examining 


the sample problem and the source listing. 


TABLE 2. FORMAT STATEMENTS FOR INPUT FILES 


1 FORMAT (11) 
7 2 FORMAT (A6,4X,8A1,2X,F19-9,19X,A5) 
3 FORMAT (3(12,3X)) 
! 4 FORMAT (13) 
5 FORMAT (8F19-3) 


aR == 


| os hata een! 


a 
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APPENDIX A 


SOME IDENTITIES IN VECTOR AND 
DYADIC ANALYSIS 


Definition A.1. A vector in R? or c? isa triple (a, 285543) 
k to denote 


> > > 
of real or complex numbers. We write A = ayi + aoi + a3 
~>s > 
that function on the set {i,j,k} of unit vectors along the x,y, and 


oe > 
z axes respectively which assigns the value a, to i, a, to j and a, 


> 


to k. 
ae : ee = 2 Gout ed ek 
| Definition A.2. If A= a,i + a,j + ak and B = by1 boJ 3 
,.! are vectors then 
>> 
AB = a,b, + aod, + ab, (A.1) 
and 
« and 
(A.2) 


| rer ® * A 
AxB = (ayb,-a 3b, )i + (azb,-a43)J + (a, b,-a5b, )k : 


= 


nae te 1 
Definition A.3. Let ww be amember of C (R?) and let A be a vector 


each component of which is a member of c! (R°). Then we define 


Tae ES 


” - UT, wey Wy (A.3) 
grad y = Vy 9x i+ ay j + 52 k 


> > 


> > > 
and for A = ai + aaj + aak ,» we define curl(A) = Vx 


Prt 
" 


aa dao > 3a) 383 + day 384 > 


}j * oe - yk (A.4) 


—-. —£)j es 


oy 9z 9Z ox 


eal 


I ROT 


Rn ac ial ee 


and 
a > ga da, da, 
div(A) Rh © ee ae ee . (A.5) 
We now introduce the concept of a dyad as follows, 
x 
Definition A.4. Suppose that the unit vectors i = (1,0,0), 
> > 


j = (0,1,0), and k = (0,0,1) are considered to be letters in an 


alphabet. Then form all the two letter "words" 


allie onaliie callie on oii en aiiieon allie on tiie on oe oo 


W= {ii,ij,ik,ji,jj,ik,ki,kj,kk} . (A.6) 


A dyad G can be thought of as a "three by three matrix" expressed in 


the form 


> — SNES eure peg 
Go = te eh tal + Fy ayid + ae es + 8¢2,1)4! + aco ayJJ Es 
rm > ae es 
S¢a,3)8% * pg yh) © Seg ahs * 8g aght + (A.7) 


Also, we may Say that G is a function from the set W into R or 1C 
at al ax als all ol oe all a oe al eo 


whose values on iji,ij,ik,ji,jj,jk, ki,kj,kk are respectively 


ad 6 PD Lind es Lila 6 Pe) ded OD an eo ee OS es Ded cal oO De 


Dyads are formed naturally in several ways. 


Settnidion Ax Me cat wee * 7 Pe 2 me on 
efiniftion A.5. If A= a,! 1 aoJ + azk and B= byl + boJ + bk 
ond a cod > a oS 
are vectors, then AB = a,b, ii + a,boij + a,bzik + aod, Ji a anboJJ + 
» ~ a >> 
+ aob,ik + a,b, ki + a,b,kj + a blkk . (A.8) 


ae, 


zo Se i 


— et ~ 


nN 


Pat ec a ae 


es 


Definition A.6. If b,,b5, an are functions in C'(R*), then 


Bs 
> > > > = 
B 


> 
B= L Vs boJ as bak implies that grad(B) = VB = 


ab, > db, > db, 
= a ae tLjlet == 


ab, +> ab, >> ab, >> 
+—— ji + ret + —— 


ab, > ab, > db, > 
Mi ered (sl eae IIR ee (A.9) 


az dz 
Whereas we used (A.5) and the divergence operation to obtain a scalar 
valued function from a vector valued function, we also use the divergence 


to obtain a vector valued function from a dyadic valued function as follows. 


> > ae ge 
Definition A.7. Let G= acy yt Tappan h eBoy rag the 9(2,1)9! 
> > Sd Rar 


F Op 9jJ) F Oyo 5,JK F Ory yy * 803 ay KJ 


a7 


+ 13, 3)** (A. 10) 


and define 


div(G) = VeG = [(a/axdary 1) + (a/ayda(o 4) + (/a2.a¢5 411 + 
+ 
[(8/ax)a(, 5) + (8/ay)a(p 9) + (a/AZ)ar3 oy] + [(3/ax)a(, 4) + 
> 
(8/8y)acz 5) + (3/8z)a(3 3y]k ‘ (A.11) 
11 


sc a he ls a lace a a el ha le a si 


Now we may make the following generalization of the Gauss divergence 
theorem which we will use later on. 


Theorem A.1. Let G be a differentiable dyadic function on a neighborhood 


of the region V of R° with boundary 9aV. 
> 
Let n denote the normal vector to the surface 3V bounding V. 


Then 


aS 
| neG da = | div(G) dv , (A212) 
aV V 


where da is the surface area and dv is the volume differential of aV 
and V_ respectively. 


Proof of Theorem A.1. The proof follows by using (A.11) and the usual 
Gauss divergence theorem component by component. 

For the sake of completeness we list the following identities for 
vector valued functions,which we also use later in the paper. 


Proposition A.1. For all differentiable vector valued functions P and 


Q we have 
> > > > > > 
div(P x Q) = Qecurl(P) - Pecurl (Q) CA. 15) 
and 
> > > 
curl(curl(P)) = grad(div(P)) - AP , (A.14) : 
> 
where the Laplacian A acts on each component of P_ by the rule, , 
52 92 92 ; 
ay = $+ 2h 428 (A.15) | 
ax dy dz 
iz 


SRD ET Os 


f 
f 
i 


E > > > 
for all scalar functions wy. Also P= b,j + boJ + b 


> 
k implies 


3 


> > > > - 
grad(P) = VP = (Vb, i + (Vb,)j + (Vb2)k » (A. 16) 


2 3 


> > > > > 


> 
div(Px cur! (Q)) = curl (P)*curl!(Q)-Pecuri(cur!(Q)), (A.17) 


> > > > 
div(P x curl! (Q) - Q x curl (P)) = 


> > > > 
Qecurl (cur! (P)) - Pecurl (curl (Q)) (A. 18) 


> > > > > > > ac 
curl (PxQ) = P div(Q) - Qdiv(P) + (Q*V)P - (P*7)Q (A.19) 


> > 
where for any vector functions A and B we define 


>> 


> > > 
(AeV)B = Lb, i te Lboj + Lbzk (A.20) 


where the operator L is defined by the rule 


Ly = a, (8/ax)p + a, (3/dy)y + a, (3/92) 


for any scalar functions jp. 


We also need the following. 
> 


Proposition A.2. Let 4 be a scalar function and A a vector 
function. Then 


> 


Q = oA 


13 


OL RCA eT Nee Eo a 


eee 


a nc hE ll A mH i a i A a sc eg cm 


5 dis: ‘oe — oo 
Sot es Bt 


implies 


VxQ = Vo x A. (A.21) 


We make use of this result and the following immediate consequence 
of the Gauss divergence theorem to derive the basic integral equation of 


electromagnetic scattering. 


Proposition A.3. Let V be a region with a boundary surface aV 


> > 
whose normal is n. Then for all smooth vector valued functions P an 
mhose normal is at for at smoormn vector valued runcrTions 

> 


Q we have the relation 


a 


| 


> > > > a 
| (Px(¥xQ) - Qx(VxP))*n da 
aV 


> > > > 
= | (Q*¥x¥xP - Pe9x¥xQ)dv . (A.22) 
V 


Remark. The reader may find it interesting to extend the definitions 


A.2 and A.3 to functions defined on the n-letter words formed from the 
>_> > 


"letters" i,j, and k. These generalizations will be the subject of 


a future report. 
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APPENDIX B 


DERIVATION OF THE FIELD EQUATIONS 
OF ELECTROMAGNETIC SCATTERING 


We use Maxwell's equations in the form 
aad” (B.1) 
(B.2) 


(B.3) 


2: D=.. (B.4) 


The constitutive relations have the form, 


> 


B= pH (B.5) 


> 


eE. (B.6) 


sources are 7) magnetic current density, i electrical current 
density, o', magnetic charge density, and pe, dectrical charge density. 
According to Stratton [6, page 229], the magnetic moment of a magnetized 
piece of matter is 


m = i] r e*(Enst)dv + fr o*(Eyn,c)da (B.7) 
V av 


> 
where w* is a surface density and r is the position vector of points in 
> 


V or aV. We define the magnetization vector M by the relation 


Ms fmdv. (B.8) 
v 


oe ss apnea aS 


SRE 


Let us assume J* = 0 andp* = 0. Also assume that B and E are time 


harmonic, which means that 
> a 
B = B(x,y,z)exp(-iwt), 


and 


E = E(x,y,z)exp(-iwt). 


Then (B.1) -(B.6) and (B.9) -(B8.10) 
imply that 


and 


V-E =0/e. 


The continuity equation implies that 
V-d - wwe = 0. 
From (B.11) and (B.12) we deduce that 
oo 2 > > 
Vx Vx E zw peE + iwypd. 
If we define 
ade t 
we deduce from (B.16) and the definition 
p= &, 
that 
= 27 > 
Vx Vx p= kE + iwyd. 
We define 
@ = exp(ikr)/r 


16 


(B.9) 


(8. 


(B. 


(B. 


(B. 


(B. 


(B. 


10) 


11) 
12) 
.13) 


14) 


15) 


16) 


19) 


en eee a 


and 

> > 

Q=¢A,; (B.20) 
where 

sg > + > 

A= A,i + AJ + A,k 


is aunit vector in an arbitrary direction. 


A calculation shows that 
A + ko = 0, (B.21) 
and 


gxyxQ = Kho +(A-v)ve . (B.22) 


Let us observe that 


> + > 
(A sv) ve =A, S(¥6)+ AotZy vo) + Aa(se V9) 


3x9 


> 
?, i [ate * anes : nie 


a*o 26 
i nes, : A2ayaz + Ag372 


= grad (Aye + Ait a Asa ) 


= v(A-vo). (B.23) 


* 2 2 2 
= ings +a,22 Asbxaa| 


Therefore, we conclude that 
> > > 
Vx 7x Q = Ak + v(A-v) . (B.24) 
Also, 


Vx Vx P= kK2E + joy: (B.25) 


17 


Thus, from proposition A-3 we see that 


f (Q-xra? - P-vx9xQ)dv = 

[t A+ (k2E + iwy 3)- E-(k2A 6 + ¥(A«v6) }dv 

- {le (x R) - oh x (oxE)] aa 

* J | @<€)-cvonty - twu( xii) -n]da 

2 {|cret ) x vol - wu (voxA)-(Ax nda 

= scot) x % + iwus (nxt) ]°A da (B.26) 


Simplifying (B.26) we find that 


[(G.vxv%P - >. vxvxQ)dv 
V 


| = Af iwydd dv - fE-v(A-v¢6) dv 
V V 


{ ~ 3 [fecre x 0+ twug(nxH)} aa (B.27) 


Now we make use of the fact that 


J(E(v-R))-n da = f div(E(v6-A) )dv 
S V 


= j[E-v(v6-a) + (¥6-A)(V-E)] dv 
V 


[ E-0(v6-A) = -A- f (v-E)vedv - f (E-n) ve aa| (B.29) 
V S 


st wnt 


ment 


Substituting (B.29) into (B.27),we deduce that 


> + > + Wa 4, 
a.) {inwd@ + (V-E)ve}dv - f (E-n)v aa] 
V S 
= Al ((nxE) «9 + fou (nH) ¢ ida] . (B.30) 
$ 
Therefore, since A is arbitrary 
> 
f Gude + (v-E) vo} dv = 


>_> > > > > > > 
J {(E-n) vw + (nxE) x V6 + iwu(nx H) 9 }da. (B.31) 
S 


Let S consist of the boundary of V together with a sphere, 9B, of radius r 
circumscribed about the point (x!,y!,z!) whose normal is directed out of V 
towards the center of the ball B of radius r,- We must compute the contri- 


bution of the right side of (B.31) to the integral over the sphere. We have 


% = r ‘kexp(ikr) P exp(ikr) @ : (B.32) 
where ep is a vector from the center of B outward. Since 
27 exp(ikr)) 
Lim ff f iw u(nxfi) r° sing do de, = 0 » (B.33) 
Mis0 00 ' : 


we observe that 


Lim 


Fi +0 


I Cin ulntte + (nx E) x 96 + (mE) ve }da 
aB 


> + > — 
= Pim | [ inxe) «06 + (n= E) Vo dda. 
3B 
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Thus, 
[tiwws + (o/e)Vo Hav = 


Y) +0 


Lim fF v(r,) | (x8) x + GER Jaw, (8.34) 
3B 


where we have used the fact that 


ory) = (b= 4k) @ (ri) 


We now prove the following Lemma. 


ia > > > > > > > >> 
Lemma B.1. If n = Ai + Bj] + Ck, then (nxE)xn+ (n-E)n = 


(A2+ B2+ C2) E 
Proof of Lemma B.1. 


Observe that 


> > > 
i(BE, ~CEY) - J(AE, - cEY) + k(AE,, - BE) 


> > > 
i (BE, ~CEY) + J(CE, - AE.) o k(AE, - BEY) ‘ 


> > 
j k 
CE,-AE. AE BE 
B C 


any 


BE,-CE 
A 


y 


Thus, 
> > > 
(nxE)xn = 


> 
2 = = 2 j 
(Cc E, ~ ACE, ABE, +B Ei 


+ (A+ Jey - BAE, - BCE,]) 
. - ((B + AYE, - CBE, - ACE, 1k : (B.39) 
We conclude that 

(nxé)xn + (me E)n = {[(B24C2)E, - A(BEy#CE,)]i 
+ [(A*+C* JE, - B(AE, +CE,)]3 
+ [(A24B2)E, - C(AE,+ BEy)]K) + [A(AE,*BE,+CE,) 13 
+ B(AE,+BE,+CE,)j 
+ C(AE, + BE, + CE, )k 

= (A24B24C2)(E 7 + EJ + Ek) : (B.40) 


This completes the proof of Lemma B.1. 


then 


— 


] 

| i Lemma B.2. If ¢(r) = exp(ikr)/r, 
? 

{ 


>> 
| Lim f (igu (nxH)o+ (nxE) x 0 + (n-E)vO}da 
i Y1 +0 3B 


= 4nE(x},yl,z}) , ro 


| r, = x-x!)2 + (y-yl)* + (z-z (B.42) 
and where B is the ball of radius r, centered at (x!,y!,z!). 


Proof of Lemma B.2. If n= ~ep = - 6/1196 Il, then applying Lemma 


B.1 and (B.33) - (B.36), we observe that 


> + > > > > + + 
Lim f {iwn(nxH)o + (nxE)x Ve + (n+ E)V¢ }da 


ke fer) E(x,y,z)da 
Pi+0 9B 


20 EE; i a 
Lim f f ik exp (fkra) exp(ikr;) E(x.y.2z)r, singdede 
0 rm 


iso r r2 
1 1 


E(x',y's2') ff sim(s)dede = 4 E(x'sy',z'), (B.43) 
0 0 


Theorem B.1. In a nonmagnetic media with : = 0 we observe that 


E(x', y', 2') = i [ (iwnde + (v-E) v6 )dv 
7 aN 


q 1 > > > + > > > > 
| -a-/ [iw u (mx tiye + (nxE) » ve + (n €) v6 |¢a : (B.44) 
S 


Proof. This follows from the basic relation (B.26) and Lemmas B.1 and 


> > > > 
B.2. where P= E and Q= A. 


{ 
| 

F | A material nonmagnetic body can be conceived as a region of free space 
| filled with charges with density ©, a current with flux density J, and 

| 


bounded by a surface with a charge density n- If we assume time-harmonic 


fields,then the continuity equation implies that 
> 
div(J) + wp = 0 (B.45) 
which implies that 
ad * -] Asan if eee 
| -( qLeiva) = daiv(d) | (8.46) 


Now we wish to define surface charge density n. In figure B.1 J- denotes 


the current density inside the surface and J+ denotes the current density 
immediately outside the surface. 


3} 


J+ 


J- 


Figure B-1. Surface charge density. 


| 
| 


4 
> + > > 


| n(ttat)A - n(t)A + [ (J, -n) A - (J_-n)AJat = 0. 


By conservation of surface charge density 


ae 


= 0 implies that 


iwn = +(J_-n) » 


or 2 
n = -i(d-n)/w > 
or 
ah fSeny 
m Tw 


(B.47) 


By the method of retarded potentials and the theory of single-layer distri- 


butions (Jones [3], p. 38, and Stratton [6], pp. 


+> 


E = rd, is 0(7) ve dv(r) 


“ahs [Jn (FY dal?) 
S 


ae le 3(F) o dv(F) . 


From (B.48), (B.47) and (B.46) we deduce that 


> 


E = 


me w 


de div(d)v¢ dv 


ev wy 


"  ¥ 


i (J:n)v da - apni J(Fe dv . 


187-188), we deduce that 


(B.48) 


(B.49) 


By dyadic analysis for all vectors A and B , we have 
+> >> > > 
div(AB) = div(A)B + A - (grad B) - (B.50) 
Thus, 


[If] div(AB)dv = ff(m-A)B ds 


V S 
= fff [aiv(ays + A-grad 3] dv. (B.51) 
V 
Letting 
A= (B.52) 
S| and 
B = grad¢ (B.53) 
in (B.51),we deduce that 
] 2 
falter Sff div(J)gradedv = 
V 
| : ] >> > 
Tne, Ri (n-J)grad.¢da - fff d-grad, graded] . (B.54) 
v 
Substituting (B.54) into (B.49) we deduce that 
é | Ex! Jf (n-d)grad,oda - fff J-grad grad odv 
| 4negiw S+r V 
q - Frcata JI (w-d)orad,oda ~ 88 17) ScHeav(r) (B.55) 
| S T y 


which implies that 


> 1 cd 
Ex « wae Ms J-grad, grad. dv 


- Ivo Sfp J(F) ¢ av(F) 
4n V 


1 


Boe) Yorn nd da, 
+ em ff (n-d)grad, ¢ 4a (B.56) 
£ 


re) 
where Ey is a spherical surface of radius p surrounding the observation 
point. 


Lemma 8.3. If o(r.r) is defined by 


o(r.r) = exP(-iko[r-r]) (B.57) 
lr-F] 
then 
sn J (n-J)grad_oda =- “E J(r) . (B.58) 
© 


Lim 5 ah 
Proof of Lemma B.3. Observe that 1+ if (n-J)grad.oda = 
z 


(2) 
a i (singcose A+ singsineB + cos¢C)grad oda, (B.59) 
a 
24 > > > 
where J = Ai + BJ] + Ck and 
be Fee = ieee ; FG 
grad @ = —5 (i singcose+ j singsine+ k cose) 


r) 


™ F(p5438)/p (B.60) 


where F(p,9¢,0) iS a bounded function. 


Thus, substituting (B.60) into (B.59), we find that 


26 


can 


Lim j; (n-J)grad,. oda = 
p+0 zr 
re) 


2m 
- fi fl, (sin ¢cos 0A + singsineB + cos gC) (singcose) singdgde 


> 27 7 


+j ff (sin cos eA + singsineB + cosgC)(singsine )singdgde 
0 
ant 

+k ff (singcoseA + singsineB + cosgC) cos¢singdgde |. (B.61) 
00 


To evaluate the right side of (B.61) we must evaluate three simple 


integrals. | 
These are 
1 2n | 
i Jf sin? ocostedgde , (8.62) | 
00 
2n 1 
1, = ee sine g¢sintedgde , (B.63) | 
0 
; and | 
2n 7 
I,= ff sin gcos* odode. (B.64) 
00 


Observe that I, = I, and that 


li et i 


T 
3 
I, = 2n(SS 4%) | = An/g - (B.65) 
Hence, 
gun { 
| Ls ae f (I-cos29)d cos ¢ { 1+ cos(26)| 4, 
> =0 0 2 
T 
3 
= ag (cos ¢ - oo $ ) . = 4n/3 % (B.66) 
: ; 
27 


—_— - 


Thus, we deduce that 


Lim f (f-S)grad, gda = -(4n/3)(Ai + 83 + Ck) (B.67) 
P+0 


Xo 


bg - 
when J is a constant vector. Since Jr) is continuous, this proves (B.58). 
Substituting (B.58) into (B.56) we obtain the following result. 


Theorem B.2 If J is the current density in a region, V, of free space, 
then, 


> 1 > 
E(r) = - Tre,To [fi J-grad grad, ¢dv 
V 
ead ff I(r) oav - cr J(r) (B.68) 
4n V ) 


at every point r in V. 


We now introduce the Green's function dyadic. 


Definition B.1. For every pair of points r and F we define 


a(t) = - 21+ na grad_grad, ) os I) (B.69) 
0 


where k* = whe and where I is the identity dyadic. 
ares iO ee ee 


: 


Theorem 8.3 If J is the current density in a region, V, of free space, 


E(r) + Hie, = juuglPv| f JCF) -G(r,#)dv (B.70) 
V 


where the right side of (8.70) is a principal value integral. 


Proof. Observe that 


iwng . ian, : 1 
: Ank? ee Amy e 4nive 
and that 
P grad. grad, @ = ~grad.grad. ¢ . 


Thus, (B.70) follows from Definition B.1 and (B.68). 


Define Gy P y(rsr’) to be the (p,q)-entry of the Green's matrix. 
p?*q 


We have 
1 32 ; 
Gy, y(ror') = 4a s + ly gee alee | (8,71) 
(Xp »Xq) 0 |°(p,q) k? aXpaXq 
Observe that if we define 
vlo) = g(ryn') = exPC-ikg|r-r'}) | (B.72) 
pear] 
and 
o = [rer'{, (B.73) 
then 
3299 
x, = v'(o) %p ~ % ) (B.74) | 
: 0 
and 
2 (x, -x2)(x,-x!) 
= — i 2 
q°"p 0” 
. | 
+ w'(p) on = ey . (B.75) . 
2 
Qe 


ea 
| 


= bs met 


Lemma B.4. 


If we assume that the penetrable dielectric is replaced 


by free space with an equivalent current distribution, then this distribution 


is given by 


Jeq = ie + jw(e(r)-e,) F- r(reE . 


Proof. Let us assume that we are in free space. 


equation becomes 
> 


> : > 
View = tee t= dy 


q° 

However, the actual equation is 
> > 

VxH-tweE = oE. 


From (B.77) and (B.78) we deduce (B.76). 


Thus, “ 
ES(r) = tr Feinty . Meath) 
ES(r) = py Jeg(r')-G(r,r')dv(r') - ie 

and 

(E = E ) = PV <(r‘) E(r')-6(r,r')dv(r') gi 


- which implies that 


1+ shu ee - PV i c(r'JE(r')*@(rgr')dv(r') = E'(r) 
V J 


0 


(B.76) 


Then the Maxwell 


(8.77) 


(B.78) 


(B.79) 


(B.80) 


(B.81) 


ae 


, 
| 
, 


Ame TERT 


APPENDIX C 


DISCRETIZATION OF THE INTEGRAL 
EQUATION OF ELECTROMAGNETIC SCATTERING 


We now carry out the discretization of the integral equation making 
full use of symmetry relations. The basic integral equation according to 


(B.69), (B.71) -(B.72), and (8.81) is 


[ 4 2, | E(r) - pv f c(r')E(r')G(ryr')dvir') = EV(r) , (c.1) 
V 

where 

G(r.r') = ~iwuo |T+(1/kg)arad, grad, |¢ (ror) : (c.2) 


Let us introduce some notation that will enable us to represent the dyadic 


in a compact form. We write 
p= |r-r'| = 7x, x1) #(x,-x5) (5-4)? , (C.3) 


exp(-ikg|r-r'|) 


vlo) = o(r.r') = (C.4) 


4n|r-r'| 
and 
’ 
oa Xi . (C.5) 
Thus, to compute grad.grad.y,we must compute (3/ax, and (3°/2x,3%4)¥ , 
Observe that 


j 
(afanp)ole) * (-tky - 5) ¥ (0) Ze 


exp(-ikgp) 
= —— vers - ax,/o” | . (C.6) 


31 


——_——<$<____—_____— 
; 


SE a Te AS 


Also appearing in the Green's function are terms of the form 


(32/ax, 3Xq)¥(0). Observe that (C.6) implies that 


a2 Ph ee (Ax, Ax 
eK = (-ik, 1/p)2( 9) | . | 


+ (1/02) (6) — + (ttre toa 


2 


- (-iky -1/)¥(p) = = 


o° 
2ik ik AX. AX 
ee Oo+ 1 + 1 + gael, pg 
| “ 0 o° p2 p2 | p* 


+ (-ikg-1/o)y(p)|(p, ‘ea - 
[-K2 + +, + et v(p) a + ae vito] Ts0) wh: (.7) 


where 5(p, q) is the Kronecker delta function. 


Substituting the definition (C.4) into (C.7) we find that 


aie « (42 + 3, + 3ikg ) exp(-ikge) axpax 
QO 


eo “ee Se 
~ (iko/p + 1/p2) CxPC-ikoo) tieah (c.8) 
To 


Another representation of (3?/ ax, 9X4) which is useful for integration is 


obtained by differentiating both sides of the expression 


mt, v'(o) (C.9) 
9 


32 


raha rps - a 


Spree 


with respect to Xq° We obtain the relation, 


AX AX 


2 ; AX, Ax 
aoa "(p) Poa “leh a re “ara | (c.10) 


Decompose the region V into subvolumes {Vi sees a Vy}. We approximate (C.1) 


by the sums, 


t(r) > N ‘ ' 
« hs Hoc, | Ely) * 2 r(r, ) E(r, ) Gri ra)V(r,) 
n¢m 


+ SM" ™ x(0)E(o)G(ryario2 singdadedo = -E"(r_) , (c.11) 


where a, is the radius of a sphere centered at Ym: 


Let us write 


Rim,n) = ht > (C.12) 
*(m, n) = k o°(m, n)? (C.13) 
and m _f 
cos of (m, ie a | (¢.14) 
rh R(myn) 
Where pe {1,2,3}. 
Then (C.7) and p = real imply that 
(1/k- 9) (87/ax p2Xq) v(e) = 
kexP tom nye 2 (m,n) (m,n) 
a ~ j > : 
“aay, msn) * 3 =(nn)*3]e08(6 x ) cos(e x; ) 
(fo(m,.n) + 085 04) . (C.15) 
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If we express (C.11) in the form 
3 N 


(myn) ¢ = 5! 
q= 2 Steger) ag" Ono vie 


then m + n implies that 


(msn)  _ -iwy,k,t(r_) 
G = 00 n = 

(V. Jexp(-i(a ) 
(xp 9X q) 4nadn n) iad (m,n) 


een 
Co 7 tees 2 Da 


+ [3 - “(m.n) + Stan) fens(ogne Yeas(aym”) ; (C.17) 


From (C.17) we deduce the relation 


gimsn) 
(x) Xq) 


~Twugk, t(r, Vi, exP(-farm ny) (a%m.n) = To(m Jn) _ 1)8(5 19) 


dna 2 
(m,n) 
. ‘ (m,n) (m,n) 
‘ 3 agin) * $0(qn)|o0s(opn n Jeostey” n)) (c.18) 
for m # n. 


We remind the reader of the relation, 


exp(-ik |r -r,|) 


o(riery) arn = (po). (C.19) 


Then, the (n,n) entry of the coefficient matrix of (C.11) has the form 
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Coane 


(n,n) _ t(r_) 
G = -§ ] alge: 
(x) 9Xq) (Pa) ; | 


: 1 a2 ' \ 
+ t(r) PV | “Tung tS g) + eee aca]? ) | dv(r'). (C.20) 
0 Pq 
n 
Here, we assume that ¥. is replaced by a sphere of radius 


/ 
i << m (c.21) 
An 


and we carry out the integration in spherical coordinates. 


We prove a series of lemmas in which we evaluate terms in the principal 


value integral defining the right side of (C.16). 


Lemma C.1. For every positive number 1 we have 


Any 27 . 
| | | 2 Le) , 2sinedgdedo = 


ny @ 6 
an 
[exe(-tk,an)-exp(-tk)] -45 | W(p)dp - (C.22) 
n 
Proof of Lemma C.1. The integral is given by 


an ft 27] , 
ae iy | usted} Zsinedgdeds = 
e 


n & 0 


arf" y'(e)odo = 4x [ eay(o), (C.23) 
n p =n 


which is the right side of (C.20). 


35 


Lemma C.2. 


For every positive number n we have 


ut A 2 ” S 
Gal X,X4(07sing)dededo 


an T 29 
mR © 6 


eee ee ee ae 
ne an 


-2 [fexp(-ikoa,) ~exp(-tkon) ] /4| 
an | 
+2 | vp do 


Proof of Lemma C.2. Observe that 


(ax;/p) = sinecose , 


(Ax%2/9) = sinesing , 


(Ax3/p) = cose - (C.25) 


Now we make use of another lemma which we use to prove Lemma C.2. 


Lemma C.3. For all p and q we have 


T 
[oxen 
aa Pz sinedgde = 65 g)4x/3 - (C.26) 


Proof of Lemma C.3. This follows from calculations used in establishing 
Lemma B.3 and from orthogonality relations. 


Continuation of Proof of Lemma C.2. The integral I representing the 


left side of (C.24) is given by 


a 
ila | “o"(o)o2de (c.27) 
n 
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Integrating by parts we find that 
a, 


n 


I= (47/3)8(, 4) pave” - | 
n 


Since by definition (C.4) and a direction calculation 


4n p p2 
r 


_ -exp(-ikge) [ ik ] 
yi(e) = bina elie ae + , 
(C.25) and Lemma C.1 imply that 
an 


2 -exp(-kgo ) 
I (41/3)5(5 q) nee cael 


-2Lexp(-ikg 9) 
4n 


an n 
+ 2 | v(e)de : 


Thus, we see that 

I = (41/3) 8¢ q) | -exe(~tkoaq) [ikoa,+1] 
+ exp(-ikgn) [ikon #11 | 
-2[exp(-ikpa,) - exp(-ikon)]/(4m) 


+2 ik val 


n 


Lemma C.4. For every positive number n we have 


a, 


ise 


(4/3) {Lexp(ikga,) - exp(-ikon)] 
4n 
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2 (it | 


[oe p*sinedededy = 


(C.28) 


(C.29) 


(C.30) 


(C.31) 


(C.32) 


Proof of Lemma C.4. This follows immediately from Lemma C.3 and 


the fact that 
an . ae 
| v'(p)pdp = ov(e)| =f v(e)de . 


n q n 


Lemma C.5. For every positive number n, we have 


I u(p)dv(o) = | anexP(-ikoan) _ exp(-ikoan) 
8(4,) ~ikg -k3 


: ait : ad 


e Ts oh as (C.33) 
Proof of Lemma C.5. We have 
ay, ay 
| v(p)p%dp = | esetcteoe lo, (C.34) 
n n 4n 
Since 
[exp(as) sds = s exp(as)/a - exp(as)/a2, (C.35) 
“n 2 1 p exp(-ikoe) exp(-ikge) | an 
| vo) odo ‘(4 jemetctten ee a ; (C.36) 
n 


n 


Theorem C.1. The principal value integral appearing as a term in the 
expression of the (n,n) entry of the coefficient matrix is defined by 


n 2 2 1 ' 
tee) = rey [-101,[5 9.4) + (1/k5) (2 /2x,9x,)] 6(tgs2 owe ) 


Vv 
n 


_ 24wuit(r,) 


= me oa © (p,q) loxP ik, a) Lika, + 1} - 1). (C.37) 
° 


Proof of Theorem C.1. We need to evaluate 


ev| (rg) tou, [6 4, + (1/2) (22/2, 344) JoCe, ") Jew ee") 


Vv 
n 


~twu,t(r, PV [ Seay + (1/2) (2/2, 344) ]oCE oe" Dav(e") (Cc. 38) 
v 


Letting BL be ¥ with a ball of radius yn removed, we use (C.10) 


to define 


Ze 
ie | coe + (1/k5) (a /2x, 2x,)] v(o)aw(e) 
B 


n 


2 ot Ax_Ax 
8 [vcrave) + (1/k°)] » (9) —P—Lav(p) 
(p,q) ° 02 


BS B 
n 


: GD 5 9) 


| YO ay (6) 


| . F Ax Ax 
-and| v0) 2g dv(p) (C.39) 
. 9 
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Using Lemma C.5, Lemma C.2, Lemma C.1, and Lemma C.4 as we proceed from 


left to right, we find that 


1 = 8(p of] exp(ik,a,) ee 


ik, -k? 
_ A n expCikyn) _ exp(-ikgn) 
| 5 


+ (1/3) 8, q) [exe(~ikg) ikon] - exptik ap) [ikgant| 
F . 1] 
- 2(exp(-ik,a,) - ooh) 


a 
+ (V8(p,q) (20/3) | Toledds 


an 


‘lee S(osa)| P(t erika ome | v (e)dp “ofa 


n 


-(1/3) [22t-hyan-er06- tka *(p oe 
0 


an 
+ (41/3) ) vole “ea (C.40) 
ny 


0 


Collecting terms in (C.40) we observe that the terms involving the singular 


integral cancel out. We have 


40 


| 
| 
| 
i 


an 


; ] in ] 
+ § exp(-ikon) |y+ - —- +--+ — 


1 
te ee (c.41) 
2 2 2 
3k, ky 3k | 


Collecting terms still further we have 


Zia 2 
| = = paslabe A eel 
n §(p,q) ase ory 3k 3k: | 


*3 e459) eo(-thon)| Si2--afe | ‘ (C.42) 


Taking the limit as n + 0 we find that 


Lim ie ge ~ Lz 
n+O - = fae |*co.a exp/( ikoa,)Likoa, +1] Ea F (C.43) 


0 


APPENDIX D 


I 
UTILIZATION OF SYMMETRY GROUPS 
IN SIMPLIFYING THE INTEGRAL EQUATION 


Define a group G generated by the operations I, Ry, Ry, and R, where 


Ry (xy »X29X3) s (-%1sX25X3)> (D.1) 
Ro(xX1.X2+X3) 2 (x1 5-X25X3)> (D.2) 
and R (x1 5%09X3) = (xy 5X25-x3). (0.3) 


Then the multiplication table for the group is given in Table 0.1. 


} TABLE N-1. THE SYMMETRY GROUP MULTIPLICATION TABLE 


I R, Ro R3 ie. Ree Ge RR 
I I Ry Ro R 3 RiR2 [| Ri Rg RoR3 Ri RoR3 
Ry Ry I RyRo Ri R3 Ro R3 Ri RoR3 RoR3 
Rp Ro ee 4 RoR3 {Ry RiRoR3 |R3 RiR3 
7 R3 R3 RiR3 RoR3 I RyRoR3 Rj Ro RyRo 
| RR, [RR OER Ry RyRoR3 {I RoR, |RiRy | Re 
| ify  [RiRs 1h RiRoR3 {Ry RoR; | I Rifts | Re 
RoRz |RoRz «|| RyRoR3 | Ry Ro Ah, jRuhe | Ri 
| RyRoRa {RyRoRz |RoRg | RRz |RiRo Rs Ro Ry I 
Let 
| G = {1,Ry,RosR3,R,Ro5R,R3sRoR3,RyRoR3}. (D.4) 


The transformation group whose multiplication group is given in 


Table 0-1 simply permutes the octants depicted in Figure D-1. 


Figure D-1. Division of space into octants. 


The basic integral equation is 


j + 2) ley - py | c(r') E(r!)-G(rye' )dv(r') = E'(r) (0.5) 


3iw Ep | : 


Assume that for all TeG, that 
G(Tr,Tr') = G(r,r') - (D.6) 


Assume that region V is invariant under a subgroup H of transformations 


inG. Then V, is the portion of V in the first octant implies that 


V = UTV, (D.7) 
TeH 
and t(Tr) = t(r) (D.8) 


for any Tc H. In what follows #(H) will denote the number of elements 


in the group H. 
Following Chen we decompose space into eight octants. Here 


Vi = {(Xsysz): x > 0, y > 0, and z > 0} is assumed to be octant one, 
Ri(Vi) is to be octant two, RiR2(Vi) is octant three, and R2(Vi) is octant 
four, octant five is R3(V1), octant six is RyR3(V\), octant seven is 
RR UR.(V),and octant eight is RR WV). Suppose that an electromagnetic 
wave is propagating in the direction of the positive z - axis, and that 
Ee! is a function of z and t only and is polarized so that E| is parallel 
to the positive x - axis. 

Suppose that 


a >; vs 


E 


n 
oO 
—s 
+ 
wn 


(D.9) 
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where et is assumed to depend only on z. Thus, 

Ch(z) = (Ez) + E'(-2))/2 (D.10) 
and 

S'(z) = (E'(z) - E¥(-z))/2 . (D.11) 
Then 

cl(z) = Ch(-z) (0.12) 
and 

si(z) = -st(-z), (0.13) 


From (0.5) and (D.7) we deduce that 


i i a E(r) - J | aH c(r')E(r')*@(ror') dv(r')I=E V(r). 
e T(V,) 


(D.14) 


Letting r' = T &', we deduce that 


=E'(r), 


; + oh len it | c(g') E(TE') G(r,Te')dv(e") 
TeH (D.15) 


Vy 


er E(sr) - J | Pv r(é') E(Ts") G(sr,Te')dv(e'=E1(Sr) - 
0 TH | 4 e 
1 m 


Now there are two possibilities. Let us define 
Be ‘ > 
C'(r) = (1/#(H)) ( F €'(sr)) (D.17) 
SeH 
and 


s!(r) = Hn C,s'(Sr), where Co = + 1 when Sr = (x,y,z) + z> 0, 
€ 


Cy = -] when Sr = (x,y,z) > z<O if r jis in the first octant| . (D.18) 


Note that if we define an operator Le by the rule 


LcP(r) = ( F(sr)/#(H) , (0.19) 
SeH 
and an operator Le in a similar way using (D.18), then 
EM(r) = Clr) + SMr) = other) + ster) (0.20) 


Lemma D.1. For every TeH, J) G(Sr,Te') = J&(Sr,e'). 
SeH SeH 


Proof of Lemma 0.1. 


Observe that 


y G(Sr,te') = J G(TTSr.TE') + 1 ere 


scH sf H (D.22) 


by (D.6). But as S runs over H, so does TS. Hence, (D.22) implies that 
(D.21) is correct. 


Theorem D.1. Let us define 


g(r) = J Eo(Sr)/#(H) 
SEH 


where E. is the solution of (D.5) obtained when F' is replaced by ¢'. Then 


| 
| h + a | v(r) - nv | ree set dv(é') = chr) . 
; 0 Se lH 


(D.23) 


Vy 


| Proof of Theorem D.1. This is an immediate consequence of Lemma D.1 


and equation (D.17) which defines (D.23). 


Lemma D.2. Let us define Cy for every T in the subgroup H of isometries 
of #% according to (D.18). Then 


L C, G(Sr.te’) = Cy F C.6(Sr,e') . (0.24) 
SeH 5 is ke 


Proof of Lemma D.2. Observe that 


») C, G(Sr,Te') = } C. G(TTSr,Te’) = J C.G(TSr,¢'). (D.25) 
SeH SeH SeH 


Since C; = 1 and Cy. = C+C,»(D.25) implies that 
| G@(Sr,Te'} #G, F Coo Gttar,c'}. (D.26) 
a: T soy 7 : 


Since as S runs over H, TS also runs over H, equation (D.26) implies 


(D.24). 
Theorem D.2. Let us define 


o(r) = | ) CE. (Sr) ih (D.27) 
SeH 


where E, is the solution of (2.5) obtained by replacing E by $'. Then 


h + FEL Joe - PV | c(e')o(E') } C.G(Sr.¢')dv or Si(r), (0.28) 
aaa SeH 
Vy 


Proof. Multiplying all terms of (D.16) by Cy and summing over S in H 
and dividing by #(H),we deduce that 


as ol t(E' ‘ TE ' ‘ 
j + she |* os Py F hermes 2 Oshlsrete ance | 
1 


: = Mr) (D.29) 
| From Lemma D.2 we deduce that 
| > 
1+ qt) jo(r) - t(é')CrE(TE') 
| + Fel he Ha PY | ‘ene 5. {o ) |dv(e') 
; 1 
= V(r) (D.30) 
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which immediately yields (D.28). 

We note now that ¥ satisfies the same integral equation that is 
satisfied by / By the uniqueness theorem we see that ‘. must be 
equal to v. Note however that the fact that (Tr) = V(r) for every 
T in H implies that the same is true of t.. This means that we can 
completely determine all values of E. just by solving an integral 
equation involving integrals over part of the first octant. 

Secondly, observe that $ satisfies the same integral equation 
that is satisfied by * By the same reasoning we conclude that 
t = $ everywhere in V. We observe that if r is in the first octant, 
then $(Tr) = c,3(r) for every transformation T in H. This 
being true of L- as well enables us to deduce the value of f. in 
any octant just by knowing its value in the first octant. 

: Finally, since E = E. + - » we see that we may determine the 
value of e in any octant simply by knowing the values of Ey and 


Ls in the first octant. 


Exercise D.1. Reconstruct the integral equation satisfied by E 
from the integral equation satisfied by ¥, 


Solution of Exercise D.1. Observe that since 5@ = [I that 
| cl(r) - { + Hel ltr) - 
i Siwe, 


-PV|{r(r* )¥(r*) ) G(Sr,r')dv(r') 
q SeH 

} ] 

| 
| 


n 


-Pv[c(r* )8(r") J G(SSr,Sr' )dv(r') 
y SeH 


1 
ae pv | <(Sr* )#(Sr" )G(r,Sr )dv(Sr') (0.31) 
SeH V 


1 


i 
| 
i 
i 


Exercise D.2. Reconstruct the integral equation satisfied by E, 
from the integral equation satisfied by $. 


Solution of Exercise D.2. We observe that since 52 = I and 
C.$(Sr') = $(r'), we see that 
Hor) = [1 + Hoefer) = 

-Pv[x(r' 8(r*) J cg6(Srer')dvir') = 
V SeH 
] 

at pv|+(Sr* )8(Sr* )og6(SSr,Sr* )dv(Sr* ) . 
€ vy, 
-Pv | (r )8(r" GU rar! )av Cr") . (D.32) 


V 
which is exactly the integral equation that was originally satisfied 
by Ee: 
Thus, Exercises D.1 and 0.2 give us proofs of the fact that 
Ee 


integral equations involving integrals over Vy the part of the scatter- 


¥ and E, = $ and shed light on the connection between the 


ing obstacle that is in the first octant, and the integral equations 


involving integrals over the entire scatterer. 


APPENDIX E 
SAMPLE PROBLEMS WITH COMPUTER RESULTS 


In this section we describe a numerical check of the program's 
validity. We approximate a ten-centimeter-diameter sphere filled 
with biological material by a complex of 528 one-centimeter cubes. 
Each of these cubes is filled with a biological material whose 
electrical properties are identical to that of the material filling 
the homogeneous sphere being approximated. We use the Mie solution 
(developed by Professor Gustav Mie in 1908) to compute the power 
density at values of r, 6, and 9 corresponding to the centers 
of the cubes (Stratton [6], pp.563-570). 


We describe here the results of two numerical experiments. The 
first experiment is scattering by a transparent region bounded by a 
sphere, and the second problem considers a translucent region 
delimited by a sphere. 


In the transparent-sphere problem a sphere with zero conductivity 
and a relative dielectric constant of one was struck by a plane wave 
with a frequency of 1000 megahertz and a strength of one volt per 
meter. The problem was solved by the computer program implementing 
the technique described in this paper. The Mie-solution method was 
also used. The y and z components of the electric vector com- 
puted by the Mie-solution program are 14 orders of magnitude 
smaller than the x component determined by this program. In the 
integral equation program described in this paper the y and z 
components of the electric vector are found to be identically zero, 
which is their correct value. The approximating cubes found in the 
first octant of three-dimensional space are called cells; there 
are 66 cubes in the first octant that are used in approximating 
the scatterer. Both the Mie-solution program and the integral 
equation program are used to compute the fields at the centers of 
these cubes; the results obtained by the two methods are listed 
in Table E-1. In this table we list the real and imaginary parts 
of the x component of the electric vector that was computed by 
the two methods. 


TABLE E-1. THE TRANSPARENT SPHERE 


INDUCED FIELD PREDICTED INDUCED FIELD PREDICTED 
BY BY 
onil INTEGRAL EQUATION METHOD MIE SOLUTION METHOD 

NUMBER Re(E,) Im(E,) Re(E.) Im(E,) 

1 9954 ~.1046 .9945 -.1046 

2 .9510 -.3093 .9510 -.3092 

3 8658 -.5004 .8658 -.5003 

| 4 7427 -. 6696 .7428 -.6695 
5 5872 -.8095 .5873 -.8094 

6 £9954 -.1046 £9945 -.1046 

‘ 7 .9510 -.3093 .9510 -.3092 
8 8658 -.5004 6858 -.5003 

9 .7428 ~.6696 .7428 -.6695 
1 10 5872 -.8095 .5873 -.8094 
it 9954 -.1046 9945 -.1046 
| 12 .9510 -. 3093 .9510 ~.3092 
13 . 8658 -.5004 .8658 -.5003 
14 .7427 -.6696 .7428 -.6695 

15 9945 -.1046 9946 -.1046 

; 16 ~9510 -.3093 .9510 -.3092 
17 8658 ~.5004 .8658 -.5003 
18 7427 ~.6696 7428 -.6695 

. 19 -9945 ~.1046 ~9945 -.1046 
| 20 9510 ~.3093 .9510 -.3092 
| 21 9945 -.1046 9945 -.1046 
22 -9510 ~. 3093 .9510 -.3092 

23 8658 -.5004 .8658 -.5003 


EE 
] 


TABLE E-1. (continued) 


. INDUCED FIELD PREDICTED INDUCED FIELD PREDICTED 
BY BY 
ent INTEGRAL EQUATION METHOD MIE SOLUTION METHOD 

veunien Re(E_) Im(E,) Re(E. ) Im(E,) 
24 7427 ~.6696 7428 -.6695 
25 5872 -.8095 873 =e 
26 9945 ~. 1046 aaa =e 
27 .9510 ~.3093 +9510 a 
28 8658 ~.5004 scat are 
29 7427. 6696 .7428 -.6695 
30 5872 -.8094 26/3 ~.8094 
31 9945 -.1046 -9945 --1046 
32 9510 -.3093 9510 -.3092 
33 . 8658 ~.5004 8658 -.5003 
| 34 7427 ~.6696 7428 -.6695 
35 9945 -.1046 9945 -.1046 
36 9510 -.3093 9510 ~.3092 
37 8658 -.5004 .8658 -.5003 
38 9945 -.1046 9945 -.1046 
39 9510 -.3093 9510 -.3092 
40 -9945 -.1046 9945 ~.1046 
41 9510 -.3093 9510 ~. 3092 
42 . 8658 ~.5004 8658 ~.5093 
43 7427 ~. 6696 .7428 -.6695 
44 9954 ~.1046 9945 -.1046 
45 9510 ~.3093 9510 -.3092 
46 . 8658 -.5004 8658 -.5003 
47 7427 -~.6696 7428 -.6695 
| 48 9510 ~.3093 9510 ~.3092 
| 49 9945 -.1046 9945 -.1046 
50 8658 -.5004 8658 -.5003 
51 9945 -~.1046 £9945 -.1046 
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TABLE E-1. (continued) 


INDUCED FIELD PREDICTED INDUCED FIELD PREDICTED 
BY BY 

eer INTEGRAL EQUATION METHOD MIE SOLUTION METHOD 

NUMBER Rett Im(E,) Re(E,) Im(E,.) 
52 -9510 -. 3093 .9510 -.3092 
53 .9945 -.1046 .9945 -.1046 
54 -9510 -.3093 -9510 -.3092 
55 . 8658 -.5004 .8658 -.5004 
56 -7427 -.6696 .7428 -.6695 
57 .9945 -.1046 .9945 -.1046 
58 ~9510 -.3093 .9510 -.3092 
59 . 8658 -.5004 .8658 ~.5003 
60 .9945 -.1046 £9945 -.1046 
61 -9510 -.3093 .9510 ~.3092 
62 .9945 -. 1046 .9945 ~.1046 
63 .9945 -.1045 .9945 ~.1046 
64 -9510 -.3093 .9510 -.3092 
65 .9945 -.1046 ~9945 ~.1046 
66 .9510 -. 3093 .9510 -.3092 


In the next problem a plane wave with a frequency of 1000 

megahertz strikes a ball whose radius is 5 centimeters, whose 
relative dielectric constant is e = 1.015625 and whose conductiv- 
ity is o = .015625 mhos per meter. The results of the comparison 
of the results of the Mie-solution program and the integral-equation | 
program described in this paper is found in Table E-2 which follows. 
In this table we list in addition to the real and imaginary parts ) 
of the x component of the electric vector the estimate of the power i 
density in each cell obtained by using the values of E and H at . 
the center of the cell. The power density is denoted by Pa in | 
Table E-2. | 
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TABLE E-2. THE TRANSLUCENT SPHERE 


INDUCED FIELD PREDICTED INDUCED FIELD PREDICTED 
BY BY 
INTEGRAL EQUATION METHOD MIE SOLUTION METHOD | 
CELL | 
NUMBER Bete Im(E,) P, Re(E,) Im(E,) a 
1 - 9376 -.0964 .00694 - 9389 -.0956 .00695 | 
2 8969 -.2941 .00690 8978 -.2943 00697 
3 .8179 -.4745 - 00698 -8181 -.4761 .00700 
4 . 7056 -.6310 -00700 . 7048 -.6336 -00702 ; 
5 -5589 -.7655 -00701 - 5639 -.7607 .00701 
6 - 9385 =-.0913 . 00695 - 9398 -.0913 .00695 
1) 8986 -.2886 - 00696 -8992 -.2901: - 00697 
8 -8205 -.4685 00697 -8201 -.4721 .-00699 
9 -7079 -.6267 00698 yoy es -.6298 00700 
10 +5596 -. 7669 -00701 - 5665 -.7574 00699 
11 - 9402 -.0814 -00695 9415 -.0830 . 00698 
12 - 9020 -.2781 - 00696 -9019 -.2819 - 00697 
13 8262 -.4562 - 00696 -8237 -.4642 . 00698 
14 wfL30 -.6163 - 00694 -7116 -.6225 - 06698 
} 15 9428 -.0673 00698 -9440 -.0704 . 00700 | 
16 - 9063 -.2658 - 00697 - 9060 -.2695 - 00698, 
17 8324 -.4468 . 00697 8291 -.4523 - 00697 | 
18 sfooe -. 5838 . 00688 - 7183 -.6113 -00695 | 
19 -9490 -.0597 - 00706 -9472 -.5345 - 00703 | 
20 9094 -.2661 .00701 Pe ig! -.2532 . 00698 | 
21 - 9380 -.0948 - 00694 - 9394 -.0924 . 00696 | 
22 8973 -.2928 - 00696 - 8988 -.2912 - 00697 
Pas, 8176 -.4739 - 00698 -8195 -.4731 00670 
24 -7071 -.6271 - 00698 . 7065 -.6309 -00701 
- 5892 -.7328 00691 -5659 -.7582 - 00700 
26 -9390 -.8994 00695 - 9402 -.0883 - 00697 


TABLE E-2. (continued) 
INDUCED FIELD PREDICTED 
BY 
INTEGRAL EQUATION METHOD 


CELL 
NUMBER 


Re(E,) 


- 8988 
-8207 
=7EE9 
- 5837 
- 9406 
- 9020 
- 8271 
«7241 
-9413 
- 9065 
- 8352 
- 9423 
-9114 
- 9390 
-8984 
- 8200 
- 6906 
-9395 
-8991 
8219 
- 7082 
+9431 
- 9032 
- 8268 
-9496 
- 9078 
«9405 
- 9001 
- 8353 
. 7048 


Im(E,) 


-.2874 
-.4671 
-.6205 
= 7379 
-.8101 
+2767 
-4533 
~.6010 
- .0638 
=.2585 
~.4277 
~.0156 
=.2261 
-.0914 
-.2893 
-.4720 
-.6510 
-.0862 
-.2857 
--4673 
—.6283 
~.0777 
=.2763 
-.4547 
-.0907 
-.2662 
-.0854 
= 52052 
-4377 
- 6364 


P 


INDUCED FIELD PREDICTED 
BY 
MIE SOLUTION METHOD 


Re(E,) Im(E,) pa 

-9001 —.2871 - 00697 
-8213 -.4692 -00699 
. 7087 = 6272 -00700 
- 5685 -.7549 - 00698 
9419 = 0799 - 00698 
- 9028 —.2789 - 00698 
-8249 -.4614 - 00698 
af32 =.6199 - 00698 
9443 -.0675 - 00700 
- 9067 =.2667 - 00698 
-8303 -.4496 - 00696 
-9473 -.0510 -00704 
-9117 = 2905 -00699 
- 9403 -.0862 - 00697 
- 9005 =. 2851 - 00698 
8221 -.4674 - 00699 
. 7098 =.6255 - 00700 
9411 -.0819 - 00698 
39019 -.2809 . 00698" 
-8239 -.4634 - 00699 
sfL2L ~.6218 - 00699 
9427 ~.0737 - 00699 
-9045 -.2728 - 00698 
8275 -.4556 - 00698 
9451 -.0614 -00701 
- 9083 -.2607 - 00698 
-9415 =.0771 - 00699 
- 9030 -.2762 - 00698 
- 8258 -.4589 - 00699 
-7146 =.6177 - 00698 


TABLE E-2. (continued) 


INTEGRAL 

CELL 
NUMBER Re(h) 

57 9403 
58 .8997 
59 .8215 
60 9409 
61 .9018 
62 9426 
63 £9423 
64 .9005 
65 £9474 

| 66 .9001 


INDUCED FIELD PREDICTED 


BY 


EQUATION METHOD 


Im(E,) 


-.0826 
=. 2891 
-.4702 
-.0409 
-.2796 
-.0704 
-.0556 
=.2905 
-.0863 
+7. IT S65 


P 


INDUCED FIELD PREDICTED 


BY 


MIE SOLUTION METHOD 


Re(E,) 


9423 
-9044 
8277 
- 9439 
- 9070 
-9461 
9429 
- 9062 
- 9438 
-9076 


Im(E,) 


~.0727 
=.2719 
~.4548 
-.0645 
—.2638 
~.0524 
-.0651 
-.2644 
-.0607 
-.2601 


A complete error analysis of the integral-equation method, not 
attempted here, will be carried out in a future report. 


In both of the numerical experiments described in this section, 
the coordinates of the outer edges of the cells is the same as that 
listed in the computer output in Appendix F, and this list is, thus, 


not repeated in this Appendix. 


Many parts of the program remain 


nearly identical to the one supplied by Professor Kun Mu Chen, and 
the authors express our gratitude to him for supplying the listing. 
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LISTING OF COMPUTER PROGRAM 


//ABEIGTGZ JOB (3H01,BC020,60,,,0010,1,¥,58) ,"HBM376@@10R COHOOR', 
// MSGLEVEL= (1,1) ,PPTY=5,TIME=56,MNSGCLASS=A 


// BXEC PLOTGCLG, REGION.GO=225K , TIME=56 
//PORT.SYSIN DD * 

REAL*2 COMP,XYZCOM, XONLYC, XANDY ,SCAT, EXPKZ,COSKZ,SINKZ,FINCX, TOTAL 
COMPLEX G(19701) ,£(198) ,C (198) , BUF (500) , FRONT, BACK, GMAT 
COMMON VOL (198,7) ,X0,¥0,20,START1,START2, STAPT3 ,ADD1,ADD2,ADD3,C1, 
102,C3,C4,FKO,NDIV,IF,10,N 

INTEGER*2 IQD,IVI,BLNK 

DIMENSICN IVI (8) ,JQ(8) ,IOD (8) 

DATA EXPKZ/"EXPKZ'/,COSKZ/"COSKZ'/, SINKZ/*SINKZ'/, BINCX/*BINCKX'/,T 
TOTAL/'TOTRL'/,BLNK/* */,XYZCOM/'X,Y,Z. '/,KONLYC/*X-ONLY */, XANDY/* X 
JANDY.'/ 

DATA IVIZ81*, 42%, 43%, 994,459,968, 078, 0887, PT73, 14159265358/,E0/8. 8 
1S4E-12/,PMU0/1.257E-6/ 

¢ READ FIPST DATA CARD FOR NDIV 
READ 5,NDIV 
5 FORMAT (I1) 
€ READ 2ND CARD FOR COMPONENTS,QUADFANTS,PREQUENCY AND TYPE OF 
c INCIDENT ELECTPIC PIELD. 
READ 10,COMP,IQD, PMEG,SCAT 
10 FORMAT (A8,2X,821,2X,P10.0,10X,A5) 

IF (COMP.EQ.XONLYC) GO TO 20 

IF (COMP.EQ.XANDY) GO TO 25 

IF (COMP.2Q0.XYZCOM) GO TO 30 

PRINT 15,COMP 

15 FORMAT ('0(',A6,') ERPOR TN COLS 1-6 OF SECOND DATA CARD') 

GO TO 395 

20 NDIM=1 
GO TO 35 

25 NDIM=2 
GO TO 35 

30 NDIM=3 

c CHECK QUADRANT CODES 

35 IP (IQD(1).EQ.IVI(1)) GO TO 45 
PRINT 40 

40 FOPMAT ('1QUADRANT CODES RRAD IN NO NOT INCLUDE QUADRANT 1°) 
STOP 

45 rsw=0 

JQ(1) =1 

DO 55 I=2,8 

IF (IQD(I).EQ.IVI(I)) GO TO 50 

JIQ(I)=0 

IP (IQD(I).NE.BLNK) ISW=1 

GO TO 55 

50 JQ(I)=I 
55 CONTINUE 
IP (ISW.EQ.9) GO TO 65 
PRINT 60,20D,30 ¢ 
60 FORMAT ('1QUADPANT CODE EERORS. SET RERD IS',10X,8(1X,A1)/24X,"SET 
1 USED IS',10X,8I2) 
65 IPPONT=JQ (5) +30 (6) +J3Q(7) +3Q(8) 

IBACK=JQ(1) #3Q(2) +30 (3) +09 (4) 

IP (SCAT.EQ.COSKZ) GO TO 75 

IF (SCAT.EQ.SINKZ) GO TO 8C 

IP (SCAT.FQ.EINCX) SO 70 85 

IP (SCAT.EQ.EXPKZ) GO TO 90 

IF (SCAT.EQ.TOTAL) GO TO 95 

PRINT 70,SCAaT 
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70 FORMAT ('0(',A6,") ERROR IN COLS 41-46 OF SECOND DATA CARD. ") 
GO TO 395 
75 ISCAT=1 
GO TO 105 
80 ISCAT=2 
GO TO 105 
85 ISCAT=3 
GO TO 105 
90 ISCAT=4 
GO TO 105 
95 IP (IPRONT.EQ.0) GO TO 100 
Iscat=5 
GO TO 105 
100 SCAT=EXPKZ 
ISCAT=4 
105 IS1=ISCAT 
IF (IS1.NE.5) GO 70 110 
ISCAT=1 
SCAT=COSKZ 
PEALD 3FD CAPD FOR NUMBEF OF CELIS IN X,Y,AND Z DIRECTIONS 
419 READ 115,NX,NY,NZ 
115 FORMAT (3 (I2,3X)) 
PEAD STH CARD FOR THE TOTAL NUMBEF OF CELLS IN THE BODY 
PEAD 120,NCELL 
120 FORMAT (13) 
MA=NDIM*NCELL 
IP (MA.LE.198) GO TO 130 
PRINT 125,MA 
125 FORMAT ('-NUMBER OF DIMENSIONS TIMFS NUMBER OF CELLS =',I&,* EXCED 
1EZS PPOGRAN DIMENSIONS‘) 
GO TO 395 
130 M1=(MA*(MA+1))/2. 
PO 135 I=1,MA 
135 E (I) =CMPLX(C. ,0.) 
PRINT TITLES 
PRINT 140 
1490 FORMAT (*1°,15X,"THE PARAMETERS OF FACH CELL AS READ IN ARE GIVEN 
IBELOW IN CNS.°/'-',6X,'N*,BX,"AMX",11X,"AMY',11X,"AMZ*,10X,'DXCM', 
110K, "DYCK", 10K, "DZCH'/) 
OMEG=2.*PI* FMEG*1. 26 
FKO=OMEG*SQRT (FMUO*E0) 
FREAD IN LOCATION AND ELECTRICAL PROPERTIES OF CELL M 
DO 170 M=1,NCELL 
PEAD 145, ANX,AMY,AMZ,RLFP1,SIG1,DKCM, DYCN, DZCH 
145 FORMAT (8F106.3) 
DYCM=DXCM 
DZCM=DXCM 
VOL (M,1) = (AMX=.5*DXCM) /190. 
VOL (M,2) = (AMY-. 5®DYCM) /100. 
VOL (M,3) = (AMZ-.5*DZCM) /100. 
DELTEX=DXCM/100. 
DELTAY=DYCH/100. 
DELTAZ=DZCM/100. 
VOL (M, 4) =DELTAX*DELTAY*DELTAZ 
VOL (M,5) =RLEP1 
VOL (M,6) =SIG1 
VOL (M,7) =DELTAX 
Z=VOL (M,3) *PKO 
GO TO (150,155,160, 165) ,ISCAT 
150 E(M) =CMPLX(-COS(Z) ,0.) 
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Go T9 170 


455 £(M) =CMPLX(0.,SIN(Z)) 117 

GO TO 170 118 

160 E(M)=CMPLX(-1.,0.) 119 

GO TO 170 120 

165 B(M) =-CEXP(CMPLX(0. ,-2Z) ) 121 

170 ‘PRINT 175,M,AMX,ANY,AMZ,DXCM,DYCM,DZCM 122 

175 FORMAT (I8,1PD15.4,5D14.4) 123 

PRINT 180 124 

180 FORMAT ('-",14X,'X,¥,Z CORRESPOND TO CENTRAL LOCATION OF EACH CELL 425 

1 IN METEPS'/'-*,6X,'N',9X,"X',13X,"Y',13X,'Z*,10X,* VOLUME", 10X, "EP 126 

1S ',8X,*STG (MHO/M) '/) 127 

DO 185 J=1,NCELL 128 

PRINT 175,3,(VOL(d,L) ,L=1,6) 129 

185 VOL(J,5) =(VOL(J,5)-1.) *OMEG*EZ0 130 

Cc PRINT THE INCIDENT FIELD VALUES AT THE CENTRAL LOCATION OP 131 

Cc EACH SMALL SUBVOLUME 132 

c PFINT PREQUENCY,MAXIMUM NUMBER OF DIMENSIONS OF BODY,TYFE OF 133 

c COMPONENTS, AND THE NUMBER OF QUADRANTS USED FOR SYMMETRY. 134 

PRINT 169,SCAT, FMEG,NX,NY,NZ,COMP,JQ,NDIV 135 

190. FORMAT ('4',9X%,* THE INCIDPNT PIELD IS OF THE FORM ',A8/*0*,17X,'F 136 

VREQUENCY =",P12.6,' MHZ.*,5X,"( NX =",13,5X,"NY =",I3,5X,°NZ =",13 137 

1,3H ) /'0',17X," FIELD COMPONENTS= ',A6,5X,*QUADRANTS USED =',812/ 138 

1°0",18X,"*NUMBER OF PAPTITIONS PER EDGE IN INTEGRATION =',I2/*0',14 139 

1X,"INCIDENT ELECTRIC FIELD HAS THE FOLLOWING COMPONENTS'/'0',12X,! 140 

| 1E INCIDENT X-DIRECTION',5X,'E INCIDENT Y-DIRECTION', 5X, *EINCIDENT 481 
1Z-DIRECTION'/*'0',7X,'N',7X, "REAL", 10X, "IMAG", 10X, "REAL", 10X, "IMAG! 142 
1,10%,*REAL*,10X,*IMAG"/) o 143 

GO TO (195,210,220) ,NDIM 108 

195 pO 200 I=1,NCELL 185 

209 PRINT 205,1,£ (I) 146 

205 FORMAT (19,2X,1P6D14.4) 147 

GO 70 230 148 

210 pO 215 J=1,NCELL 149 

JPN=J+NCELL 150 

215 PRINT 205,3,2(J),E(JPN) 161 

GO TO 230 152 

j 220 DO 225 K=1,NCELL 153 
a KPN=K+NCELL 154 
: KNN=KPN4NCELL 155 
225 PRINT 205,K,E(K),E(KPN) ,# (KNN) 156 

230 GRID=NDIV 157 
GFACT=0.5*((1./GRID)-1.) 158 

C1=-2, *OMEG*FMU0/ (2. *PKO*FKO) 159 

C2=3./(4.*PI) 160 

C3=3.*OMEG*EO 161 

C4=-OMEG*FMUO*PKO/ (4, *PI*GRID** 3) 162 

Cc CALCULATE THE G MATRIX 163 

31=0 1€8 

Ti=1 165 

DO 270 IP=1,NDIM 166 
DO 270 M=1,NCELL 167 

XO=VOL (M,1) 168 
i YO=VOL (M, 2) 169 
| Z0=VOL (M, 3) 170 
K=(IP-1) *NCELL+M 171 

DO 270 TO=IP,NDIM 172 

N1=1 173 
; IP (IP.JEQ.TO)NI=N 176 


DO 27C N=N1,NCELL 
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L=(IQ-1) ®NCELL+N 
DELTAX=VOL (N, 7) 
DELTAY=DELTAX 
DELTAZ=DELTAX 
STAP?T 1=GPACT*DELTAX 
STAPT2=GPACT*DELTAY 
START3=GPACT*DELTAZ 
ADD1=DELTAX/GRID 
ADD2=DFLTAY/GFID 
ADD3=DELTAZ/GRID 
GO TO (235,240,245) ,I0 
235 BACK=GMAT (JQ(1)) +GMAT (JQ (2)) +GMAT (JQ (3) ) +GMAT (JQ (4) ) 
PRONT=GMAT (JO (5) ) +GMAT (JQ (6) ) +GMAT (JQ (7) ) +GMAT (JO (8) ) 
GO TO 250 
280 BACK=GMAT (JQ(1)) -GMAT (JQ (2)) +GMAT (JQ (3) ) -GMAT (JQ (4) ) 
PRONT=GMAT (JQ (5) ) -GMAT (JQ (6) ) +GMAT (JQ (7) ) -GMAT (JQ (8) ) 
GO TO 250 
245 BACK=GMNAT (JQ(1)) -GMAT (JQ (2) ) ~GMAT (JQ (3) ) GMAT (JQ(4)) 
FRONT=-GMAT (JQ (5) ) +GMAT (JQ (6) ) +GMAT (JQ (7) ) -GMAT (JQ (8) ) 
259 GO TO (258,260,255,255,265) ,IS1 
255 G(I1)=BACK+PRONT 
GO TO 270 
260 G (i1) =BACK-PRONT 
GO TO 270 
265 G(I1)=BACK+FPONT 
3123141 
BUF (J1) =BACK-FRONT 
IF (J1,LT.500) GO TO 270 
WRITE (°) (BUP(I) ,I=1,500) 
J1=0 
270 Ii=I1+1 
IP (IS1.NE.5) GO TO 280 
IP (J1.FQ.0) GO TO 275 
WRITE (9) (BUP(I) ,7=1,01) 
275 REWIND 9 
CALCULATE THE INDUCED FIELDS 
280 CALL MATPCH(G,E,MA,1,TEB, 1.E-35) 
IP (IER.NE.0) GO TO 395 
PRINT 285 
285 PORMAT ('-',28X," (INDUCED ELECTRIC FIELDS IN QUADRANT 1)'/ '0',16X 
1,°N*,9X,*EX (V/8) *,11X, "BY (V/M) *,11X,"EZ(V/M) °,5X,"PWR DEN (W/CU. BM 
4FT.)*/) 
CALL PRSUB(E,NDIM,NCELL,POWT) 
CHECK POF CEEATING TOTAL POWER 
IP (IS1.£0.5) GO TO 305 
IF (IFRONT.NE.C) GO TO 395 
IF (IBACK.EQ.10) GO TO 290 
IP (IBACK.EQ.1.0R.IBACK.GT.5) GO TO 395 
XNUL=2. 
GO TO 295 
290 XMUL=4. 
295 POWT=XMUL*POWT 
PRINT 300,POWT 
300 FORMAT (*-TOTAL POWER IN THE BODY =',1PE15.5," WATTS") 
GO TO 395 
305 SCAT=SINKZ 
DO 315 I=1,MA 
C (I) =F (I) 
IF (I.GT.NCELL) GO TO 310 
E (I) =CMPLX (0. ,STN (VOL (1,3) *FKO) ) 
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310 
315 


329 
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360 
365 
370 


375 


380 


385 
390 


395 
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GO TO 315 
E (I) =CMPLX (0. ,0.) 
CONTINUE 


PRINT 190,SCAT,FMEG,NX,NY,NZ,COMP,JQ,NDIV 

GO TO (320,330,340) ,NDIM 

DO 325 I=1,NCELL 

PRINT 205,I,£(TI) 

GO TO 350 

DO 335 J=1,NCELL 

JPN=J+NCELL 

PRINT 205,J,2(3),E(JPN) 

GO TO 350 

DO 345 K=1,NCELL 

KPN=K+*NCFLL 

KNN=KPN+NCELL 

PRINT 205,K,E(K) -E(KPN) ,z (KNN) 

I1=I1-1 

K2=9 

K1=K2¢1 

K 2=K2+500 

IF (K2.GT.I1)K2=21 

READ (9) (G(K) ,K=K1,K2) 

IP (K2.1T.11) GO TO 355 

CALL MATPCH (G,E,MA,1,IER,1.E-35) 

IF (IER.NE.Q) GO TO 395 

PRINT 285 

CALL PRSUB(E,NDIM,NCELL, POWT) 

DO 360 I=1,MA 

G (I) =E (1) +C (I) 

E (I) =2 (1) -C (1) 

PRINT 365, (JQ(I) ,I=5, 8) 

FORMAT (*1FIELDS IN FRONT LAYER (OCTANTS',UI2,* )") 
PRINT 370, FMEG,NX,NY,NZ,COMP,JQ,NDIV 

FORMAT (1HO,24X,* (INDUCED ELRCTRIC FIFLDS IN ONE OCTANT) '/'O',17X, 
1" PREQUENCY =",F12.6,' MHZ.",5X,'( NX =",I3,5X,"NY =",13,5X,"NZ =', 
123,3H ) /*0',17X," FIELD COMPONENTS= ',A6,5X,"QUADRANTS USED =',8I 
12/°0',18X,* NUMBER OF PARTITIONS PEF EDGE IN INTZGRATION =",12/'-', 
116X,'N*,9X, "EX (V/M) *,11X, "EY (V/M) *,11X,"EZ(V/M)',5X,"PWER DEN (W/CU 
1. MET.) '/) 

CALL PRSUB(E,NDIM,NCELL, POWPF) 

PRINT 375, (JQ(I),I=1,4) 

FORMAT (‘1FIELDS IN BACK LAYER (OCTANTS',4I2,* )*) 

PRINT 370,FPMEG,NX,NY,NZ,COMP,JQ,NDIV 

CALL PRSUB(G,NDIS,NCELL, POWB) 

IF (IPRONT+IBACK.EQ.3€) GO TO 385 

IPF (IFRONT.EQ.11.AND.IBACK.£Q.3) GO TO 380 

IP (IFFONT.FQ.12.AND.IBACK.EQ.4) GO TO 380 

IP (IPRONT.EQ.13.AND.JQ (5) .EQ.5.AND.IBACK.EQ.5) GO TO 380 
IP (IFRONT.NE.5.OR.IBACK.NF.1) GO TO 395 

XMUL=1. 

GO TO 390 

XMUL=2. 

GO TO 390 

XMUL=4. 

POWT=XMUL*® ( POWF+ POWB) 

PRINT 300,PO"T 

STOP 

END oe 

SUBROUTINE PRSUBEL,NDIN, NCELL,PTOT) 

COMPLEX X(1) 
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COMMON VOL (198,7) 
PTOT=0.D¢ 
DO 20 M=1,NCELL 
EPA=CABS (X (M) ) 
IF (NDIM.EQ.1) GO TO 5 
J=M4NCELL 
EPB=CABS (X (J) ) 
IF (NDIM.FQ.2) GO TO 10 
K=J+NCELL 
PPC=CABS (X (K)) 
GO To 15 
EPB=0.D0 
FFC=0.D0 
EWE=C.5*VCL (M,6) * (EFA*EPA+EPB*2FB+EPC*EPC) 
ETOT=PTOT+PWR*VOL (M,4) 
PRINT 25,M,EPFA,EFB, EFC, PWR 
PORMA™ (I18,1P4D18.5) 
PRINT 30,PT07 
FORMAT (1H-,18X,*TOTAL POWER ABSORBPD IN THE PIFST OCTANT =', 1PE14 
1.5, WATTS'/'- N',GX, *EXREAL',6X,*EXIMAG', 6%, *EX-ABS', 3X, "THETA- 
1X',3X,EYREAL',6X,*EYIMAG',6X,"EY-2BS',3X,"THETA-Y', 3X, 'BZREAL* , 6X 
1, ‘EZIMAG',6X,*EBZ-ABS', 2X, *THETA-Z'"/) 
PRINT REAL, IBAGINARY,ABS. VALUE AND PHASE ANGLE OF INDUCED 
PIELD 
DO 70 M=1,NCELL 
GO 70 (45,40,35) ,NDIM 
K=M+2*NCELL 
FPC=CABS (X (K)) 
CALL ANGLES (XK (K) , EF 3) 
J=M4+NCPLL 
EPB=CABS (X (J) ) 
CALL ANGLES (X(J) ,EF2) 
EPA=CABS (X (M) ) 
CALL ANGLES (X(M) ,2F1) 
GO TO (50,60,65) ,NDIM 
PRINT 55,M,X(M) ,EPA,2P1 
FOPMAT (I%,3(1P3E12.4,CPF7.1)) 
GO TO 70 
PRINT 55,M,X(M) -2FA,EP1,X(J) ,BFB,EF2 
GO TO 70 
PRINT 55,M,X(M),E2A,EF1,X (J) ,EFB,FF2,X(K) ,BFC,EF3 
CONTINUE : 
FETURN 
END 
FUNCTION GMAT (ZIQD) 
DIMENSION U (3) 
COMPLEX TAUN, FCTR1, FCTR2;EXP0 ,TOTAL, GMAT 
COMMON VOL(198,7) ,X0,Y¥0,20,START1,START2, START3,RDD1,ADD2,ADD3,C1, 
1€2,C3,C4,FKO,NDIV,ID,10,N 
IP (IQD.NE.0) GO TO 19 
GMAT=CMPLK (0. ,0) 
GO 70 8C 
GO TO (15,20, 25,30, 35,40,45,50) ,IQD 
U (1) =X0-VOL (N, 1) 
U (2) =YO-VOL (N, 2) 
U (3) =Z0-VOL (N, 3) 
GO TO 55 
U (1) =X6+VOL (N,1) 
U (2) =¥0-VOL (N,2) 
U (3) =Z0-VCL (N, 3) 


THIS PAGE IS BEST QUALITY PRACTICABLE 
FROM OOPY FURAISHED TODDC ___- 


GO TO 55 17 

25 (1) =XC+VOL (N,1) 18 

U (2) =Y9+VOL (N,2) 19 

U (3) =Z0-VOL (N,3) 20 

GO TO £5 21 

30 0 (1) =xO-VOL(N,1) 2 

0 (2) =YO+*VOL (N, 2) 23 

0 (3) =Z0-VOL (N, 3) 24 

GO TO 55 2 

35 0 (1) =x0-VOL (N, 1) 26 

1 (2) =YO-VCL (N, 2) 27 

0 (3) =Z0+VOL (N,3) 28 

GO To 55 29 

40 0 (1) =xX0+VOL (N,1) 30 

0 (2) =YO-VOL (N,2) 31 

U (3) =Z0+VOL (N, 3) 32 

GO To §5 33 

BS 011) =x¥9+VOL (N,1) 34 

0 (2) =¥0+VOL (N, 2) 35 

1 (3) =Z0+VOL (N, 3) 36 

GO To 55 37 

50 U (1) =XO-VCL (N, 1) 38 

U (2) =YO+VOL (N, 2) 39 

U (3) =Z0+VCL (N, 3) 40 

55 VN=VOL (N,4) 41 

TAUN=CMPLX (VOL (N,6) , VOL (N,5)) a2 

IP (0 (1) *0 (1) +0(2) *0 (2) +U (3) *U (3) .NB.0.) GO TO 60 43 

IP (IP.NE.IQ) GO TOS.- 44 

FCTR1=CMPLX (0.,C1) *TAUN 45 

A= (C2*VN) **, 33333 33333333333 46 

EXPO=CEXP (CMPLX(0.,-FKO*A) ) 47 
FCTR2=EXPC*CMPLX(1.,FKO*A) -CMPLX(1.,0.) 48 

CMAT=PCTE1*FCTP 2-CMPLX(1.,C.) ~TAUN/CMPLX (0. ,C3) 49 

l GO TO 80 50 
60 TOTAL=CMELX (0.,0.) 51 

UIMIN=U (1) +START1 52 

U2MIN=U (2) #+START2 53 

1 (3) =U (32) #+START3 54 

} DO 75 IZ=1,NDIV 55 
U (1) =U1KIN 56 

DO 70 IxX=1,NDIV 57 

U (2) =U2MIN 58 

DO 6£ IY=1,NDIV 59 

B=SQRT (U (1) *0 (1) +U (2) #U (2) +0 (3) *U (3) ) 60 

A=PKO*R 61 

BSQ=A*®A 62 

FCTRI=CEXP (CMPLX(O0.,-A)) /(ASQ*A) 63 

PCTR2= (0 (IP) *U (10) / (F*P) ) *CMPLX (3.-ASQ,3.¥*A) 64 

IP (IP. 2Q.1Q) FCTR2=FCTR2-CMPLX(1.-ASQ,A) 65 
TOTAL=TOTAL+PCTP1*FCTR2 66 

‘ 65 0(2)=C0 (2) +ADD2 67 
7€ U(1)=0(1)+ADD1 68 

| 75 0 (3) =U (3) #ADD3 69 
GMAT=CHPLX(0.,CU*®VN) *TAON*®TOTAL 70 

| 80 RETURN 71 
r END 72 
SUBROUTINE MATPCH (A,B,N,M,IEP,EP) 1 

COMPLEX A(1),B(N,M),S 2 

TER=0 3 

WM1=N-1 4 
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P4=4 5 
T2=N 6 

DO 20 I=1,NM14 7 
C=CABS (A(T1)) 8 
IF (C.LT.1.0-15) PRINT 5,I,A (31) 9 

S FORMAT ('OFOR ROW',T4,* PtvoT IS*,1P2D12. 4) 10 
IF (C.LT.2P) GO TO 25 11 
IPt=I+1 12 
T3=I1 13 
K1=T34N-2414 1& 

DO 15 J=IP1,N 15 
I3=13+1 16 
S=A (13) /A(T1) 17 
J1=13 18 

DO 1 K=0,N 19 
A(K1)=8 (K1) -S*A (91) 20 
J1=314+4 21 

19 K1=K1+1 22 
DO 15 L=1,M 23 

1S B(J,L) =B(J,L) -S*B (1,1) 24 
I1=I1+I2 25 

20 I2=r2-1 26 
C=CABS (A (K1~1)) 27 
IPF (C.LT.1.D-15) PRINT 5,N,4 (K1-1) 28 
T=N 29 
IP (C.GE.EP) GO TO 35 30 

25 IER=1 34 
PRINT 30,N,I 32 

30 FORMAT ('-THIS SYSTEM OF ORDER*,I4,' IS SINGULAP FOP ROW',T4,*. 978 33 
1 THE PROGRAM CANNOT HANDLE THIS CASE.') 34 
FETORN 5 

35 I3=(N* (N41) ) 72 36 
I4=2 37 
DO 50 I=1,N 38 
K=N-1+1 39 
KP1=K+1 4&0 

DO 50 L=1,m 44 
$=0.D0 42 
IP (K.EQ.N) GO TO 45 43 
J1=T3+41 44 

DO 40 J=KP1,N 65 
S=S+A (J1) *B(3,L) 6 

4O 31=3141 47 
45 B(K,L)=(B(K,L)-S) /A (13) 46 
I3=I3-14 4s 

50 I¥=T441 50 
RETURN 51 
END 52 
SUBROUTINE ANGLES (X,EF) 1 
COMPLEX X 2 

c DETERMINE ANGLY BETWEEN REAL AND IMAGINARY PARTS OF FFIELD 3 
A=REAL (X) 4 
B=AIMAG (X) 5 

TE 4B) 6S, 45,25 6 

5 IPF (A.NE.0.) GO TO 10 7 
PBF=-99.0 8 

GO TO 30 3 

10 PR=57.2957795131*ATAN2 (B,A) 10 
GO TO 30 4 

15 IP (A.GE.C.) GO TO 20 12 
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EF=-180. 13 
GO TO 30 14 
20 FF=0. 45 
GO TO 30 16 
25 IPF (A.NE.O.) GO TO 19 17 
BF=90. 18 
30 RETURN 19 
END 20 


//GO. PTO9PCO1 DD UNIT=SYSDA,SPACE=(TRK, (10,5)) ,DISP= (NEW, PASS) 
7/GO.SYSIN DD * 


1 
X,Y,2. 12345678 1000. TOTAL 
S 5 5 
66 

1. 1. 1. 1.915625 0.015625 1. 1. 1. 

1. bs 2. 1.015625 6.015625 1. 1. 1. 

1. 1. 3. 1.015625 6.015625 1. 1. 1. 

1. 1. 4. 1.015625 0.015625 1. 1. 1. 

1. 1. 5. 1.015625 0.015625 1. 1. 1. 

1. 2. 1. 1.015625 0.015625 1. 1. 1. 

1. 2. 2. 1.015625 0.015625 1. 1. 1. 

1. 2. 3. 1.015625 0.015625 1. 1. 1. 

1. 26 4. 1.915625 0.015625 1. 1. 1. 

1. 2. 5. 1.015625 6.015625 1. 1. 1. 

1. 3. 1. 1.015625 0.015625 1. 1. 1. 

1. 3. 2. 1.015625 0.015625 ie 1. 1. 

1. 3. 3. 1.915625 0.015625 1. 1. 1. 

1. 3. 4, 1.915625 0.015625 1. 1. 1. 

1. 4, 1. 1.015625 0.015625 1. 1. 1. 

1. 4. 2. 1.915625 0.015625 1. 1. 1. 

1. 4, 3. 1.015625 ¢€.015625 1. 1. 1. 

1, 4. 4. 1.015625 6.01£€25 1. 1. 1. 

1. 5. 1. 1.015625 ©,.015625 1. 1. 1. 

1. S. 2. 1.015625 0.015625 1. 1. 1. 

2. 1. 1. 1.915625 0.615625 1. 1. 1. 

2. 1. 2. 1.015625 0.015625 1. 1. 1. 

2. 1. 3. 1.015625 6.015625 1. 1. 1. 

2. 1. 4. 1.015625 0.015625 1. 1. 1. 

2. 1. 5. 1.995625 6.615625 1. 1. 1. 

{ 2. 2 1. 1.915625 0.615625 1. 1. 1. 
2. 2. 2. 1.015625 6.615625 1. 1. 1. 

2. 26 3. 1.015625 0.015625 1. 1. 1. 

2. 2. 4. 1.015625 ¢.015625 1. 1. i. 

2. 2. 5. 1.015625 9.015625 1. 1. 1. 

2. 3. 1. 12015625 0.015625 1. 1. 1. 

2. 3. 2. 1.01562€& 9.015625 1. 1. 1. 

2. 3. 3. 1.018625 0.015625 1. 1. 1. 

2. 3. G, 1.015€25 0.015625 1. 1. 1. 

2. 4. 1. 1.015625 0.015625 1. 1. 1. 

2. 4. 2e 1.015625 ¢.C15625 1. 1. 1. 
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